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1 .   I ntroduct ion 

Divide  and  conquer  (DC)  methods  have  been  developed  for  "the  symmetric 
e igenprobl em ,  see  Cuppen  [Cu]  ,  Dongarra  and  Sorenson  [DS]  ,  and 
Krishnakumar  and  Morf  [KM] ,  and  have  for  these  problems  been  shown  to  be 
efficient  on  parallel  computers  and  competitive  on  single  processor 
machines  [DS]  ,  [Cu]  ,  [KM]  .   The  DC  method  has  also  been  applied  successfully 
to  the  computation  of  singular  values  by  Jessup  and  Sorensen  [JS] .   In  the 
present  paper  we  describe  a  DC  method  for  the  unitary  e igenprobl em ,  and  we 
also  discuss  the  simplifications  that  arise  for  real  orthogonal  matrices. 

Let  H  6  CnXn  be  unitary.   Then  H  is  unitarily  similar  to  a  upper 
Hessenberg  matrix  with  real-valued  non-negative  subdiagonal  elements.    If 
a  subdiagonal  element  vanishes,  then  the  eigenproblem  splits  into 
e igenprobl ems  for  smaller  upper  Hessenberg  matrices.   We  therefore  may 
assume  that  H  is  a  upper  Hessenberg  matrix  with  positive  subdiagonal 
elements.   Then  all  eigenvalues  of  H  are  simple.    It  is  easily  seen  that  H 
can  be  written  as  a  product  of  n  Givens  reflectors  G:  €  CnXn, 

II  =  H(7i,72..--.7n):=  GXG2  .  .  .  G^Gn  ,  (1.1) 

where,  for  1  <  k  <  n, 


Gk  :  = 


k-1 


^k   7k 


n-k-1 


>  7k  €  C  ,  <7k  £  R  ,  crk  >  0  , 

'k  I 


l7k|2  +  *k2  =  1, 


(1 .2a) 


and 


Gn 


n-1 


•7n 


?  7n  €  C,   |  7n  |  =  1 


(1.2b) 


Here  L  denotes  the  j  x  j  identity  matrix.   The  7=,  1  <  j  <  n,  are  the  so- 
called  Schur  parameters  of  II,  and  7=  denotes  the  complex  conjugate  of  1  =  . 
The  cr-  ,  1  <  j  <  n,  are  said  to  be  complementary  parameters  of  H,  and  are 
the  subdiagonal  elements  of  H. 

The  DC  method  described  uses  the  product  representation  (1.1)  of  II, 
the  so-called  Schur  parametric  form  of  H.   An  application  of  particular 
interest  to  us  is  the  computation  of  Pisarenko  frequency  estimates  for  a 
random  stationary  stochastic  process,  see  below.   In  this  application  H  is 
defined  by  its  Schur  parametric  form.   The  determination  of  Gaussian 
quadrature  rules  on  the  unit  circle,  so-called  Gauss-Szego  quarature 
rules,  also  gives  rise  to  unitary  (or  real  orthogonal)  matrices  in  Schur 
parametric  form,  see  Section  5.   When  the  Schur  parameters  are  not 
explicitly  known,  they  can  be  computed  from 


Ti  =  "H 


11 


tj  =  -cg; 


1  Gj-2 


G2  G"  H)jj'   J  =  2'3<  •  •  •  >n 


where  G.   denotes  the  conjugate  transpose  of  Gk  ,  and  M ^  denotes  the  element 
(j,j)  of  a  matrix  M  G  Cnxn  . 

The  DC  method  is  most  easily  described  for  II  G  RnXn  orthogonal .   Then 


H  =  G1    G: 


Gs.!  Gs  Gs+1 


Gn 


II 


In- 


G« 


1, 


(1-3) 


where  Hj  G  RsXs  and  H2  €  R(n"s)x(n"s)  are  orthogonal.   The  Givens  reflector 


Gs  G  R    ,  s  <  n,  can  be  written  as  a  Householder  transformation 


G,  =  I  -  1\ 


(1.4) 


wh 


ere 


w 


Wc 


es+1ws+i     6 


1/2 


(1     +     78  ) 


1/2 


;s+l 


:=    -2-1/2(l     -    7s)l/2 


(1.5) 

(1 .6a) 
(1.6b) 


Throughout  this  paper  e:  denotes  the  jth  column  of  an  identity  matrix  of 
appropriate  order.   By  (1.3)-(1.4),  H  is  orthogonal  ly  similar  to 


H, 


(I  -  2wwH)  =:  H  -  2HwwH 


(1.7) 


This  is  one  step  of  the  DC  method  for  orthogonal  matrices:   The  eigen- 
values, and  if  so  desired  the  eigenvectors,  of  H^  and  H2  are  computed 
first.   H   is  a  rank-one  modification  of  H,  and  the  eigenvalues  of  H  are 
computed  as  the  zeros  on  the  unit  circle  of  a  rational  function,  whose 
poles  are  the  eigenvalues  of  H  . 

Section  2  describes  the  DC  method  for  unitary  matrices.    In  Section  3 
we  show  some  results  on  the  orthogonality  of  the  eigenvectors  and  on  the 
location  of  the  eigenvalues.   These  results  are  analogous  to  bounds 
presented  by  Dongarra  and  Sorensen  [DS]  for  the  DC  method  for  symmetric 
matrices.   Section  4  discusses  simplifications  that  can  be  made  when 
H  is  real  and  orthogonal,  and  also  considers  some  computational  details. 
Computed  examples  for  the  orthogonal  eigenvalue  problem  are  presented  in 
Sect  ion  5 . 

An  outline  of  a  unitary  DC  method  with  convergence  results  for  the 
root-finder  has  previously  been  presented  in  [GR]  .   The  splitting  into 
subproblems  is  done  differently  in  the  present  paper.   A  related  DC  method 
is  described  by  Arbenz  and  Golub  [AG]  .   Cybenko  [Cy]  reduces  the 
orthogonal  eigenproblem  to  an  eigenproblem  for  a  symmetric  tridiagonal 


matrix.   The  orthogonal  eigenproblem  is  in  [AGR1]  solved  by  solving 
singular  value  problems  tor  certain  bidiagonal  matrices,  and  a  CJK 
algorithm  tor  unitary  matrices  is  presented  in  [Grl] .   A  comparison  with 
respect  to  accuracy  and  speed  of  these  methods  still  remains  to  be  done. 
Here  we  only  note  that  DC  methods  are  suitable  tor  implementation  on 
parallel  computers,  see  [DS] , [KM] ,  and  Section  2.   Some  of  the  schemes 
mentioned  transform  the  orthogonal  eigenvalue  problem  to  a  symmetric  one. 
The  real  eigenvalues  of  the  latter  problem  are  then  mapped  to  the  unit 
circle  to  yield  the  eigenvalues  of  the  orthogonal  eigenproblem.   The 
mapping  from  the  interval  to  the  unit  circle  may  be  sensitive  to 
perturbations  of  arguments  near  the  end  points  of  the  interval ,  and  it  may 
therefore  be  difficult  to  determine  eigenvalues  close  to  ±1  accurately 
with  these  schemes. 

Pisarenko  [Pi]  proposed  a  method  for  decomposing  a  random  stationary 
stochastic  process  {xm}??  q  ,  xm  G  IR  ,  into  a  sum  of  harmonics  in  white 
no  i  se  ,  i.e., 


^m  =  Z  Qp     cos(m^  +  ^)  +  ym,  m  >  0,  (1-8) 

£=1 


where  the  6„    are  arbitrary  phase  shifts  and  {ym}°°_n  is  a  zero  mean  white 
noise  process  with  variance  a    .   The  </>„    are  called  P  i  sarenko  frequency 


't 


est  imates .   Assume  for  simplicity  that  p  is  the  number  of  distinct 
harmonics  in  the  'signal'  {xm}°°_Q  is  known,  and  that  0  <  </>,,  <  n    for 
1  5:  ^  <  P-   Then  the  <j>»    can  be  determined  as  follows.   Form  the 

(2p+l )  x  (2p+l )  Toeplitz  covariance  matrix  M  for  the  signal  {xm}°"?  q  ,  and 

2  2p 

compute  its  least  eigenvalue  Amjn .   Then  Amjn  =  a    ,  see  [Pi]  •   Let  {->,}•  i  be 

the  Schur  parameters  associated  with  the  Toeplitz  matrix  M-A  f  I.   They  can 


be  determined  from  "the  Szego  recursions  (Levinson's  algorithm) ,  see  e.g 
[AGR2]  .   From  our  assumptions  it  follows  that  M-Amjr)l  is  singular,  but 
leading  principal  submatrices  not  identical  with  M-AmjnI  are  not. 


Therefore,  -1  <  7j  <  1  for  1  <  j  <  2p ,  and  72p  <E  {-1,1}.   By  (1.1)-(1.2)  it 

'ollows  that  the  Schur  parameters  {7:} 
.2PX2P  ...^.u  j;„4.!«-*  „; „!..-„  /\\2P 


follows  that  the  Schur  parameters  {7:}-_ i  define  an  orthogonal  matrix  H  € 


^z\>~*v  wj-th  distinct  eigenvalues  {A:}.  ,  .  Enumerate  the  eigenvalues  so  that 
those  with  Im(A:)  >  0  have  smaller  index  than  the  eigenvalues  with  Im(Aj)  < 
0.   Then  the  Pisarenko  frequency  estimates  are  given  by 

</>j  :=  arg(Aj)   ,     1  <  j  <  p  . 

The  coefficients  Q:  of  (1.8)  are  two  times  the  weights  belonging  to  the 
Gauss-Szego  quadrature  rule  with  abscissas  A:,  1  <  j  <  p.   For  details  see 
[AGR2] ,  where  also  references  to  related  work  can  be  found.   The  unitary 
DC  method  yields  the  Gauss-Szego  weights  with  no  extra  computational 
effort  when  computing  the  eigenvalues  Aj  .   Gauss-Szego  quadrature  is 
discussed  in  Example  5.1  of  Section  5. 


2.   The  unitary  e igenprobl em 

In  this  section  we  describe  a  divide  and  conquer  method  for  unitary 
right  Hessenberg  matrices  with  positive  subdiagonal  elements.   First  we 
need  to  generalize  the  splitting  (1.3)-(1.7)  to  Givens  reflectors  with 
complex-valued  Schur  parameters.   This  is  accomplished  by  noting  that  the 
Gs  defined  by  (1.2a)  are  diagonally  unitarily  equivalent  with  a  real 
Householder  transformation.   Introduce 


7s  := 


7s/ I  7s  I   ,      7s  #  0 


7s  =  0 


Then     |  7s  I     —    1?     and    for    Gs    defined    by     (1.2a)     we    obtain 


Is+1 


7s 


G« 


il 


s-l 


n-s-1 


hs  I       ^s 
^s        I  7s  I 


n-s-1 


(2.1) 


where,     similarly    to     (1.5)-(1.6),     w    =    esu;s    +    es  ,  jws  ,  1     G    Rn    and 


cs     :=    2-l/2(l    +     |7s|)l/2 


u;s+1     :=    -2"1/2(1     -     |7s|) 


1/2 


(2.2a) 
(2.2b) 


Substitution  of  (2.1)  into  (1.1)  yields 


H  = 


II. 


In-« 


(I  -  2wwH) 


(2.3) 


ith 


Hi     :=    H(7j  ,72,  .  .  .  ,7s-i  >  "7s')  , 


■=-     '-  zr     I 


-     I 


H2     :=    H(7S  7s+i  »  7s  7s+2  '  •  •  •  '7s  7n) 


A  unitary  similarity  transform  of  (2.3)  yields,  analogously  with  (l.r), 


H'  :  = 


Hi 


H, 


(I  -  2wwH)  =:  H  -  2HwwH  . 


(2.4) 


Let 


Hk  =  Wk  Ak  WR!  ,   k  =  1,  2, 


(2.5) 


be  spectral  resolutions,  i.e.,  the  Wk  are  unitary  and  the  Ak  are  diagonal. 
Then  H  has  the  spectral  resolution  H  =  W  A  W  ,  where 


W  :  = 


W, 


,  A  := 


=  diag(Ax  ,A2 ,  .  .  .  , An)  , 


(2.6) 


with  |  Ak  |  =  1  for  1  <  k  <  n. 

We  are  in  a  position  to  describe  how  the  spectrum  of  H  can  be  obtained 
from  A,  the  last  row  of  Wj  and  the  first  row  of  W2 •    Introduce  the 
characteristic  polynomial 

X(A)  :=    det(AI-H)  =  det(AI-H')  =  det (AI -H+2HwwH) 
=  det(AI-H)  det(I+2(AI-H)"1HwwH) 
=  V»(A)(1  +  2wH(AI-H)"1fiw) 


.Hr, 


1aiT/H. 


=  0(A)  (1  +  2wnW(AI-A)"1AW"w)  , 


8 


where  H'  is  defined  by  (2.4),  W  and  A  by  (2.6),  w  by  (2.2)  and  xl>(\) 
det(A[-H] ) .   Let  z  =     [G\=i    be  given  by 


=  WHw  = 


W^es^s 
W2elws+1 


(2.7) 


and  define  the  spectral  function 


<j>(\) 


*(A)  _  1      ,  o„H/XT   n-U„  _  1   ,  o  £  \,   |2 


^(A) 


+  2zM(AI-A)-1Az  =  1  +  2  £   |  <•  | 

j=l       A  "  AJ 


n      „  A  +  A: 

j=i      A   AJ 


(2.8) 


u 

/here  we  have  used  that  z  z  =  1 .   Let 


9j  :=  arg(Aj)   ,   6     :=  arg(A)   ,   0  <  0j  ,  0  <  2* 


Then,  with  i  :=  4  - 1 , 


n       _      ,Q,    -    6, 
<KA)  =  i  E   ICjI2  cot(-i-2— )  =:  i*(0)  , 


$ 


,<?=  -  0, 


'CO    =   i.E     iCjlV-inf^J!)  >    i^   =   i 


(2.9) 


(2.10) 


We  may  assume  that  the  6     are  distinct  and  that  all  (,-     ^  0,  because 
otherwise  we  can  make  these  conditions  hold  by  deflation,  see  below.   Let 
6-     G  [0 ,  2?r  [ ,  1  <  j  <  n,  denote  the  zeros  of  $(6).       Then  the  eigenvalues  of 
H'  and  of  H  are  given  by  Aj'  :=    exp(i0j'),  1  <  j  <  n.   The  sets  {^}H=1  and 
{#:  }:_i  strictly  interlace. 

We  describe  a  rootfinder  for  $>($).       By  the  inequality  (2.10).  the 
zeros  of  $(#)  can  be  determined  accurately.   We  may  assume  that 


0  <  01    <    02    <     •  •  •  <  #n  <  27r  and  that  9        ,  our  initial  approximation  of  a 
zero  of  $(#)  5  satisfies  #n-27r  <  0  <    0-^  .        By  the  strict  interlacing  of  the 

sets  {#:}P_i  and  {0-   }"— i  >  ^(^)  has  precisely  one  zero,  denoted  0-^    ,  in  the 
open  interval  ]0n-27r,  0±  [ .   Assume  for  the  moment  that 

$(0(o))  <  0  ,  (2.11) 

and  introduce 

*(0)  :  =  p  +  a    cot^1  2  -)  .  (2.12) 

The  coefficients  p  and  a    are  determined  by  osculatory  interpolation,  i.e., 
$(0(o))  =  $(0(o))  ,   l>'(0(o))  =  $'(^(o))  ,  (2.13) 

wh i  ch  y  i  e Ids 

t,    =  *</o))  -  *'(#(0))  sin(«!  -  0(o))  , 

I  =  af <#«•>)  sinf  i  -/°'))   .  I 

The  zero  0    of  $(#)  in  ]^n-27r,^1[  is  our  next  approximation  of  0-^   .       New 
approximations  0  of  0j  are  computed  from  0      ,  j  >  1,  in  a  similar 

fashion.   The  sequence  {0      }^n  satisfies  0         <    0^     for  j  >  0,  and  converges 
monoton i cal ly  and  quadrat i cal ly  to  0±     as  j  increases,  see  [GR]  for  a 
proof . 

If  instead  of  (2.11)  we  have 

<J>(0(O))  >  0  ,  (2.14) 

then  we  replace  (2.12)  by 

*(0)  =  p    +  a    cot(gn  ~    g)  (2.15) 
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in  (2.13).   This  defines  p    and  a .   The  zero  8  of  $(6)     in  the  open 


l  nterva 


1  ^On-lit  ,6±  [  is  our  next  approximation  of  0^'.        New  approximation: 


,(j+l) 


i(J) 


of  #j   are  computed  from  0    for  j  >  1  in  a  similar  fashion.    II 


le 


sequence  {#  }\2.n    satisfies  9         >    6±     for  j  >  0,  and  converges  monoton i cal 1 y 
and  quadrat i cal ly  to  6^   ,  see  [GR] .   In  the  implementation  used  to  generate 
the  computed  examples  of  Section  5,  the  iterations  are  carried  out  until 
<J>(0(j+l))  >  $(0(i)  <  ^{6{}'l}).       The  value  0(j)  is  accepted  as  an  approximate 
root  of  $(0)  =  0. 


From 


6   ' 
A  :  =  d  i  ag  [e    ,  e  ',...,  e  n  ] 


W'       Wr' 


(2.16) 


and  the  spectral  resolutions  (2.5)  of  H x  and  H2 ,  we  can  now  compute  the 
spectral  resolution  of  H: 


H  =  WAW 


H 


WMW  =  I 


(2.17) 


By  (2.3) -(2.4) ,  we  can  for  some  vector  u  6  Cn  express  H  as 


H  =  H  -  2uuH  , 


u  :  = 


H^es^s 


elu's+l 


e  cn 


Let  A  :=    exp(i#)  and  v  6  Cn  form  an  eigenpair  of  H,  i.e.  Hv  —    vA .   Then 


.-^H 


(H  -  2uun)v  =  vA  , 


or,  equ i val ent ly , 


.~.H. 


(H  -  IA)v  =  ua  ,     a  :=  2unv  . 

This  shows  that  any  normalized  eigenvector  v  of  H  associated  with  the 
eigenvalue  A  is  a  normal ization  of 


1  1 


V 


=  (H-AI)_1u  = 


-i, 


(H1-IA)-xH1esws 


-l 


(H2-IA)-ie1wi+x 

W1(I-AI1IA)-1WI1,esu;s 
W2(A2-IA)-1\\fesu;s+1 


(2.18) 


where  Wk  and  Ak  are  given  by  (2.5)  .  Let  ||  ||2  denote  the  Euclidean  vector 
and  matrix  norms.  From  HWjj  =  ||W2||2  =  HAJ^  =  1  and  (2.6),  (2.7),  (2.10), 
(2.18)  it  follows  that 

6(A)  :=  ||v'||2  =  H(A-IA)-1  x2  ||2 


n 
=   •  1      -   I2 


)1/2  =  (J*'(*))l/2  >  J 


(2.19) 


We  choose 


VA  = 


W1(l-AI11A)'1w5Iesa;s 
W2(A2-IA)-1W^e1cs+1 


/6(A) 


(2.20) 


and  note  that  the  lower  bound  (2.19)  for  6(A)  indicates  that  severe 
cancellation  of  significant  digits  does  not  take  place  in  the  computation 
of  vA  by  (2.20)  . 

By  (2.7) ,  we  only  require  the  last  row  of  Wj  and  the  first  row  of  W2 
(as  well  as  A)  in  order  to  determine  the  spectral  function  (2.8) .   Hence, 
if  we  do  not  desire  the  eigenvectors,  then  it  suffices  to  determine  the 

first  and  last  rows  of  W  in  order  to  be  able  to  compute  the  spectrum  of 

H     H 
larger  problems.   We  call  the  triplet  {A,ej W,enW}  the  partial  spectral 

resol ut  ion  of  H.   The  first  and  last  elements  of  v,  can  easily  be 

determined  from 


e?vA    =    e[IW1(I-A}IA)-1w5Iesu;s/6(A)      , 
■e|?vA    =    en.sW2(A2-IA)-1W^e1o)s+1/6(A) 


(2.21) 
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We  may  assume  that  the  columns  o"f  W  are  such  that  all  components  of  the 


vector  W  et  are  real  and  positive.   Then  e.  W  e1     is  the  weight  correspond- 
ing to  the  node  exp(i#k  )   in  the  Gauss-Szego  quadrature  rule  with  nodes 
{exp  (  i  6-')  }"_i  5  see  [Gr2]  and  Example  5.1. 

We  assumed  above  all  components  (,*    of  z  to  be  non-vanishing  and  all 
eigenvalues  A:  of  H  to  be  distinct.   These  conditions  can  be  made  to  hold 
true  by  def lat  ion .   Our  discussion  follows  Dongarra  and  Sorensen  [DS]  . 
First  assume  that  Q  vanishes.   By  (2.4)-(2.7) , 

WHH'W  =  A  -  2WHHwwHW  =  A (I  -  2zzH)  ,  (2.22) 


and  from  e*  z  =  0  it  follows  that 


A(I  -  2zz  )  e  p    —    Ae«  =  A<,e«,     \„     :=  e„  Ae 


I    /vc£ 


(2.23) 


Substituting     (2.23)     into     (2.22)     yields 


H'We„    =    We.A 


rt 


and    therefore 


W, 


W2A2' 


eC    = 


W, 


W2A2r 


eCXC 


Thus     if    £»    =    0,     then    we    can    determine    an    eigenpair    of    li    without    explicitly 
computing    a    zero    of    $(8)     and    without    using     (2.20)  . 
For    a    general     z    €    C    ,     with    z    z    =    1 ,     we    obtain 

||A(I-2zzH)e£||2    =    2|C£|      , 
and    we    accept    {A^,e/>}     as    an    eigenpair    of    A     if 
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2|C£l<<i  (2.24) 

for  some  smal 1  constant  (j  . 

Assume  that  z  =  [Ci]"_i  with  2  |  C;  I  >  (i  for  all  j.   We  may  now  be  able  to 
deflate  due  to  close  eigenvalues.   Let  A  —    diag[Al5A2,  .  .  •  ,  An]  with 
A,  «  A2  ,  and  choose  the  Givens  reflector 


G  = 


so  that 


-7  a 
a        7 


n-2 


€  C 


nxn 


(2.25) 


Gz   =    [(Kit2    +    IC2I2)1 /2,o,C3,C4-  .  .  .  ,Cn]T    , 


1  .  e 


'»    :=    IC2l/(Kil2  +    IC2I2)1/2    , 

■7    :=    -Cx    ^/(ICil2   +    IC2|2)1/2 


2     ^H 


We    accept    {7^      +    72  I  7  I     » G    e2}    as    an     (approximate)     eigenpair    of 
A(I     -    2zzH)     if 

|7(r(A1     -    A2)  I     <    t2, 


(2.26) 


for  some  smal 1  constant  e2 ,  because 


M\  oH 


2npH. 


A(I-2zZn)G"e2    -     (X^    +    72l7r)G"e2||2    =     |7<r(7l    -    72)  I 


If    7j    =    72    then    we    have    determined    an    eigenpair    exactly.        In    case 
Aj    7^    A2    we    note    that     |  y<r  ( Aj    -    A2)|     <    i|Aj    -    A2|,     and  ,     moreover ,     if     |  7  |     «    0 
or     I7I     »    1    then     ^^(Aj    -    A2)  |     <<     |  X1    -    A2  |  .        Hence,     inequality     (2.26)    may 
be    satisfied,     even     if     |  Xl    -    A2  |     >    e2  •       Assume    that     (2.26)     is    valid.       Then    A 
is    replaced    by 
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A     :=    diag[A1(72)     +    \2a2  ,  A3  ,  A4  ,  .  .  .  ,  An]     G    C<n-1>*("-1) 

and  if  A  has  close  eigenvalues,  "then  deflation  is  repeated. 

The  unitary  DC  method  can  be  used  in  two  ways.   One  approach  is  to 
divide  the  original  e igenproblem ,  as  well  as  subproblems  so  obtained, 
until  only  trivial  e igenprobl ems  of  orders  two  and  one  remain.    These 
small  e igenprobl ems  are  solved  analytically.   From  the  solutions  of  small 
e igenproblems ,  the  solutions  of  e igenprobl ems  of  larger  size  are  computed, 
and  this  step  is  repeated  until  the  solution  of  the  original  eigenproblem 
has  been  determined.   This  approach  is  used  in  the  numerical  examples  of 
Sect  ion  5 . 

An  alternative  approach  is  to  use  the  DC  technique  to  generate  just  a 
few  sube igenprobl ems ,  each  of  which  can  be  solved  independently  by  some 
other  numerical  scheme,  such  as  the  unitary  GR  method  [Grl] ,  or  the  scheme 
in  [AGR1] ,  in  case  the  matrix  is  real  orthogonal . 

We  conclude  this  section  with  some  bounds  of  the  computational 
complexity  of  the  unitary  DC  method.   Assume  that  H  £  CnXn  is  given  in 
Schur  parametr  i  c  form  (1.1)  with  pos  i t  i ve  subd  i  agonal  elements  a-  .     Let  n,= 
2   for  some  positive  integer  C,  and  subdivide  the  given  eigenproblem  unl i 1 
5  e igenprobl ems  for  2x2  matrics  are  obtained.   The  latter  e i genprobl ems 
are  solved  analytically.   We  assume  that  the  number  of  iterations  required 
by  the  rootfinder  for  $(#)  can  be  bounded  independently  of  n. 

Let  first  n  independent  processors  be  available.   The  reduction  of  the 
original  eigenproblem  for  H  to  ^    e i genprobl ems  for  2x2  matrices  can  be 
carried  out  in  tj  :=  0(log2S)  time  steps.   This  computation  only  requires 
the  determination  of  the  Schur  parameters  for  the  unitary  matrices  of  the 
smaller  e i genprobl ems ,  see  (2.3).   Let  the  Schur  parameters  for  all  ^ 
unitary  2x2  matrices  be  known.   The  spectral  resolution  of  all  these 
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matrices  can  be  computed  in  t2  :  =  0(1)  "time  steps.   Assume  that  the 
partial  spectral  resolutions  (2.21)  of  all  2      unitary  2^   x  2^   matrices 
are  known  for  some  j  €  [2 , £]  .    In  order  to  compute  the  partial  spectral 
resolution  of  all  2    unitary  2/  x  z  matrices,  we  have  to  compute  2/  zeros 
of  each  of  the  2    functions  $(#)  ,  see  (2.8)  .   Hence,  a  total  number  of  n 
zeros  have  to  be  computed,  and  we  use  one  processor  to  determine  each  one. 
Each  function  $(#)  has  2/  terms,  and  can  therefore  be  evaluated  in  0(2/) 
time  steps  for  each  value  of  6.        Hence,  we  can  determine  all  eigenvalues 
of  all  2    unitary  2^  x  2^  matrices  in  to   :=  0(2;  time  steps.   For  each 
eigenvalue  we  compute  the  first  and  last  elements  of  the  corresponding 
eigenvector  from  (2.21).   The  first  and  last  element  of  one  eigenvector 
can  be  determined  by  one  processor  in  0(2^)  time  steps.   These 
computations  have  to  be  carried  out  for  n  eigenvectors  by  n  processors  and 
therefore  require  t^   =  0(2^)  time  steps.   Hence,  the  number  of  time  steps 
required  to  determine  the  partial  spectral  resolution  of  H  by  n  processors 
i  s 


*1  +  t2  +  £  UJ;  +  £  Vr    =    0(n) 
j=2        j=2 


(2.27) 


Now  let  n   independent  processors  be  available,  and  assume  that  the 
partial  spectral  resolutions  of  all  2      unitary  2J"1  x  2J"1  matrices  are 
known  for  some  j  £  [2  ,  C]  .   We  have  now  n  processors  available  for  each 
evaluation  of  each  of  the  2    functions  $(0).        Each  of  these  functions  $(0) 
has  2r  terms  and  can  for  each  value  of  6    be  evaluated  in  0(log22^)  time 
steps.   Hence,  we  can  compute  all  eigenvalues  of  all  2    unitary  "2:   x  T 
matrices  in  to  i—    0(log22^)  time  steps.   The  first  and  last  elements  of 
each  eigenvector  of  each  of  the  2  J  unitary  2*  x  21  matrices  can  be 
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determined  in  t^   :=  0  (  1  og2  2^)  time  steps,  by  using  n  processors  to 
compute  each  sum  with  "2r     <    n  terms.   The  initial  determination  of  the  R 
unitary  2x2  matrices  and  their  spectral  resolutions  cannot  be  sped  up 
essentially  by  using  more  than  n  processors,  and  requires  {^(X*^)     +    0(t2) 
time  steps.   Hence,  the  number  of  time  steps  required  in  order  to  compute 
the  partial  spectral  resolution  of  H  by  n   processors  is 

0(ti)  +  0(t2)  +  £  t!.j)  +  £  tj,j)  =  O(log2n)  +  £  0(j)  =  0(log2n)  _   (2.28) 

j=2       j=2  j=2 

The  time  complexities  (2 . 27) - (2 . 28)  suggest  that  the  unitary  DC  method 
presented  could  be  attractive  for  use  in  real-time  signal  processing 
appl i  cat  ions. 
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3.   Some  properties  of  the  unitary  DC  method 

We  show  some  properties  of  the  eigenvectors  of  II  and  zeros  of  $(#)  . 
Analogous  results  have  previously  been  obtained  by  Dongarra  and  Sorensen 
[DS ,  Lemmas  4.2,  4.6  and  4.7]  for  the  DC  method  for  the  eigenproblem  for 
symmetric  tridiagonal  matrices.   The  following  formulas  are  used  in 
several  of  the  proofs: 


*(A)  =  1  +  2  £      Kj|2  t—^x  , 
j=l       A  "  Aj 


i  n         o        A: 

*'(A)  =  "2  £   Kjl2 J— 

j=l       (A  -  Aj 


(3.1) 


(3.2) 


$'(0)  =  <^'(A)A 


(3.3) 


Lemma    3 .  1       Let    Aj/x    G    C,      |A|     =     | /i  |     =1    and    A    ^    n .       Assume    that    Aj/i    £    {A:}",  , 
where    A    =    diag[A1?A2,  .  .  .  ,  An]  .        Let    v,     and    v^    be    defined    by     (2.20)  .       Then 


|vA%|     =     U'(A)^(/i)|-l/2 


4>{\)     -    <j>(n) 


=    |*'(<?A)*'(^) 


■1/2 


A    -    ft 


(3.4) 


e  -    e 


A         lCV 
where  A  =  e   ,  //  =  e    ,  0  <  9^  ,  6  ^    <    2n     .    In  particular,  if  A  and  /x  are 

distinct  eigenvalues  of  H,  then  <f>(\)    =    <f>(n)     =  0,  and  therefore  vAv„  =  0. 


Proof.   By  (2.20),  (2.6)  and  (2.7), 


VA  = 


W, 


W. 


In-s 


(A  -  IA)"1Z/6(A)  , 


(3.5) 


and  therefore 
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'a%  -  siha     nylCA     ^yl  *oo 


=     (SiX^SOOY1!:      n Ufr r     •  (3.6) 

j=i    (Aj    -    A)  (Aj   -   /<) 


Now 


AAj  x     /    A;  A; 


(Aj-A)(Aj-/i)  (Aj-A)(Aj-,0 


hbfi  ~  v-J  •  '<3-7) 


Substituting  (3.7)  into  (3.6)  yields 


i    a  ^  #•    KjI'aj     o  £    LiLjV 


vA"v,  =  (a.(A)»W)-»  ^(2  t      ^    2±      ^)  , 


and  by  (2.19),  (3.1)  and  (3.3), 


vA%  =  (^(A)KaOA,0-i/2a  ^(A^  :  lM 


,              ,                i/o    i6\    *(0\)  -  $(^u) 
=  ($'(^)$'(^))-l/2  ie  A  -^ ^ 


e    -  e 
This  shows  (3.4) .     D 

The  denominator  |A-^|   in  (3.4)  suggests  that  it  may  be  numerical ly 
difficult  to  obtain  orthogonal  eigenvectors  when  the  associated 
eigenvalues  are  very  close.   The  following  lemma  sheds  some  light  on  this 
situation,  and  shows  that  due  to  deflation  the  roots  of  $(#)  are  bounded 
away  from  each  other. 

Lemma  3 . 2   Let  Ai  =  e   ,  1  <  j  <  n,  be  the  eigenvalues  of  A,  and  let  z  = 
[Cj]jl_i  be  defined  by  (2.7).   Let  (l     be  an  arbitrary  but  fixed  positive 
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constant  and  assume  that  the  Aj  are  pairwise  distinct  and  that  2  |  G  |  >  Cj 
for  all  j.   These  conditions  can  be  made  valid  by  deflation.   Assume  that 
the  A:  are  sorted  so  that  0  <  0X    <    02    <     •  •  •  <  0n    <  27r  ,  and  let  0n+1  :  =  lit    + 
61.        Let  0    G  [0,27r[  be  a  zero  of  $(9),    and  let  k  be  such  that  0k    <    0    <    0k  +  1. 
Then 


,(3.8) 
(3.9) 


Proof.    Introduce  the  index  sets 


Ix     :=    { j  :     0    <    0,    +    2tt£    <    0    +    n ,     for    some    I    G    Z,     1     <    j     <    n}     , 
I2     :=    {j  :     0-tt    <    0j    +    2;r£    <    0,     for    some    £    G    Z,     1     <    j     <    n}     . 


Then     Ix     D     I2    =    0    and     Ix    U    I2    =    {1 , 2 , . . . , n} .        Further 


A  +  A; 
1  A 


+  \-  0    -    0        (   <    c 

=->j  =  «*HH)  |  >  0 


0    0:x  f  <  0  ,      j  G  Ij  , 

j  e  i2 


In  particular,  k  G  I2  and,  provided  that  k  <  n,  k+1  G  Ix.    If  k  =  n  then 
1  G  Ij.   Moreover, 


,0-0 


-      0: 


:ot( 2"^)  ^  cot(— 2— j)  ,        Vj  G  Ij 

(0    -    0kx        ,0    -  0;n 
:ot[ 2— 5J  >  cot( 2— Jj  ,  Vj  G  I: 


(3.10) 


From    $(0)    =    0,     it    follows    that 


,0-0; 


E     Kjl2  cot(V)  =  £     |Cj|»  cot(V) 
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(3.11) 


By  (3.10),   (3.11)  and  z  z  =  1  we  obtain,  provided  that  k  <  n, 


-KR+1|2  cotf4*±i)  <  c„t(H)  , 


(3.12) 


or,  equ i val ent ly , 


Kk  +  i|2  tan(^)  <  tan(  "  +  ,1   ) 


(3.13) 


If  k  =  n  then  we  define  Cn+l 
Now  assume  that 


(,'n  and  (3  .  1 2)  -  (3  .  13)  remain  valid. 


(3.14) 


We  wish  to  determine  a  lower  bound  for  #k+l-^-   From  tan  ^  <  x  for 
0  <  x  <  ^  it  follows  that  if  0  <  0k+1-#  <  %f ,  then 

tan(gk+12"  0)  <  ek+1  -  e   . 

Substituting  (3 . 14) - (3 . 15)  into  (3.13)  yields 

|Ck  +  1|2  tan(I(flk  +  1  -  *k))  <  <?k  +  1  -  6     , 
and  from  tan  (^  (^k  +  l~^k )  )  -  aC^k-fl-^)'  we  obtain 

Kk  +  ll  4(^  +  1"  °k)     ^  ^k  +  l  "  e     ■ 

Finally,  substituting  |<k  +  i  |  >  ^  into  (3.16)  yields  (3.8). 

In  order  to  show  (3.9),  we  note  that  from  (3 . 10) - (3 . 1 1 )  and  z  z 
f ol 1 ows  that 

-  e, 


(3.15) 


(3.16) 


=    1    it 


-    c 


o*f-^±!)    >    |Ckl>    cot(^) 


or,     equ i val ent ly , 


tan 


re  -  e 


0..J-,   -  *■ 


f-^)>    |Ck|2tan(!^__!)   , 


(3.17) 
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which  corresponds  "to  (3.13)  .   We  now  assume  "that, 

*k  +  1  -  9    >    ±(0k  +  1  -  0k)     .  (3.18) 

We  would  like  to  determine  a  lower  bound  for  #-#k  .   Similarly  as  in  the 
derivation  of  (3.15)  we  obtain  that  if  0  <  #-#k  <    =j2r  ,  then 

2~JS)  -  9    "  ^  •  (3.19) 

From  (3.  17)-(3.  19)  and  tan  (J  (0k  +  1-0))  >  l(0k  +  1-0)  ,  we  obtain 

'  -  *k  >  lCkl24(^k+l  "  *k)  •  (3'20) 

Finally,  substituting  |  Ck  |  >  U     into  (3.20)  yields  (3.9).      D 

Our  final  lemma  shows  that  the  computed  eigenvectors  are  close  to 
orthogonal  if  the  zeros  of  $(#)  are  evaluated  with  sufficient  accuracy. 

Lemma  3.3.   Let  A  =  diag[A1,A2,  .  .  .  ,  An]  ,  and  let  A,  /2  be  computed 
approximations  of  the  distinct  roots  A,  /j    of  <f> .    Introduce  the  relative 
errors  ak  ,  /?k  of  Ak-A  and  Ak-/i,  respectively,  i.e. 

rAk  -  A  =    (Ak  -  A)(l  +  ok) 

<  k  =  1,2, . . . ,n  .  (3.21) 

Uk  -  A  =  (Ak  -  /i)(l  +  /?k)  • 

Assume  that  for  some  constant  0  <  e  <  1,   |  o-k  |  <    (    and   |  /?k  |   <  e  for  all  k, 
and  that  |A|  =  |  ft  \     -    1.   Then 


ivi%i  =  iVcv^i  <  «(2  +  o(Hi)2  • 


where    C    =    diag[/7l5p2?  •  •  •  j /?n]     with 
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Pk  '■ 


ak    +    /?k    +    Qk/?k        6(\)6(fl) 


(1    +    ak)  (1    +    /3k)\6(A)6(/x) 


(3.22) 


For    rj     :=    e       ,  0  <  On    <    2tt,     we    define    6  (?;) 


=  (^'(r/)^)l/2  =  (J*'(^)) 


1/2 


Proof.  We  first  note  that  since  A,  p.  are  computed  by  determining  zeros  of 
$(0)>  0  <  6  <  2n  ,  the  requirement  |  A|  =  | /i  |  =  1  i  s  satisfied.  Analogously 
to  (3.6)  we  obtain 


v.%  =    (S(X)SWr1   t     - — N 

A  "  k=l  (Ak  -  A)(Ak  -  A) 


-l 


(6(\)6(ii)y'  e 


Clcl 


k=i    (*k  -  ~X)(K  -  /i)(i  +  sk)(i  +  /?k) 


where  the  last  equality  follows  from  (3.21).   Now 
0  =  VA/J  =    (^(A)6(m))'1  £ 


Ki 


k=l    (Ak  -  A)  (Ak  -  p) 


and     (2.19)     imply    that 


Kkl2 


n 

k=i    (Ak  -  A)  (Ak  -  *0 


=    0 


Therefore 


A       /^  \,.  =  1  (Ak-A)  (Ak-/i)  (1+Qk)  (l+/?k) 


ickl2 


n 

k5i  (\-    v 


K 


=  (l(X)'(wrl£  (wvrt((^)W)     !)  ' 


(3.23) 
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which  shows  that  v-Hv-  =  v^Cv^  with  C  defined  by  (3.22).   From  (2.19), 
(3.3),  (3.2)  and  (3.21)  it  follows  that 


k=l 


Ki 


2  (-Ak)A 
(A-Ak)2 


(3.24) 


n   / 

E  (id 

k=i  v 


2   ("\c)A       A/A 


(A-Ak)2   (l+ok) 


From  (-AkA)/(A-Ak)2  >  0  it  follows  that  A/(A(l+ak)2)  >  0  and  therefore 
VA     ^  ,„   ,   ,_  IN-2  ^  ^   ,  ^-2 


(1  +  «k) 


>   (1  +   lOkl)"^  >  (1  +  0 


(3.25) 


Substituting  (3.25)  into  (3.24)  yields  <5(A)/i5(A)  <  1  +  c,  and  similarly  one 
can  show  that  6(ji)/<5(/i)  <  1  +  c.   Hence,  by  (3.22)  , 


Pk\     1 


(      +      €      +      €' 

(i  -  o2 


(1  +  e)2  =    e(2  +  0(r^l): 


Fi  nal ly , 


lvIHvZil     =     lvAhCv^l     ^    HvaII2    HCII2    IM2    =    IIC||2    =   max        |  pk  | 
^       ^  "  1 <k<n 

and  the  desired  bound  now  follows  from  (3.26) .      D 


(3.26) 
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4.   The  orthogonal  eigenproblem 

The  computational  work  required  for  the  real  orthogonal  eigenproblem 
is  smaller  than  for  the  unitary  one.   Th i s  section  discusses  these 
differences,  and  considers  some  details  of  our  implementation  of  a  DC 
scheme  for  the  real  orthogonal  eigenproblem.   Our  computer  program  is  for 
the  case  when  H  G  Rnxn  with  n  =  2  ,  where  £  is  a  positive  integer,  and  we 
assume  in  this  section  that  n  is  of  this  form.   Many  of  our  comments  are 
valid  for  more  general  values  of  n,  also. 


We  first  note  that  the  subdivision  of  the  eigenproblem  for  II  into 
smaller  e igenprobl ems ,  as  described  by  (1.3)-(1.7),  does  not  require  any 
computational  work.   Subdivision  yields  the  block-diagonal  matrix 


H  :_  G1G3G5  .  .  .Gp.sGn^Gn^Gn  , 


(4.1) 


and  we  obtain  simple  formulas  for  the  eigenpairs  of  each  2x2  block  on  the 
diagonal  as  follows.   Let 


G  :  = 


-7  a 
a        7 


►2X2 


-1  <  7  <  1,  <r    >    0,  72  +  a2    =    1 


(4.2) 


Since  G  is  real,  symmetric,  orthogonal  and  has  distinct  eigenvalues 
{^1*^2}  >  we  have  Aj^  =  1  and  \2    —    -1.    Let  Xj  =  [(p^]    be  an  eigenvector  of 
unit  length.   Then  we  can  choose 


^  =  .2-l/2(l  +  7)"1/2 
X2    =  2"1/2(1  +  7)1/2, 


(4.3) 


ind    from    <r    =    (1    -    72)1         it    follows    that 


i,    =    2-l/2(l  7)1/2 

,*2    =    ,2"1/2(1  7)"1/2 


(4.4) 
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Cancellation  of  significant  digits  is  avoided  by  using  (4.3)  if  7  >  0  and 
(4.4)  otherwise.   An  eigenvector  associated  with  A2  =  -1  is  given  by  x2  :  = 

[£2>-£l]T- 

If  7n  =  -1  then  Gn  =  ln,     and  the  eigenpairs  of  Gp.jGp  are  those  of  Gn_±  . 

If  7„  =  1  we  need  to  determine  the  eigenpairs  of 


-Tn-l  -o-n-1 
^n-l   "Tn-l 


We  find  that  the  eigenvalue  Aj  =  _7n-l  +  i^n-l  °^  ^  has  an  associated 
eigenvector  xx  :=  [2"    ,  -2"   ]   ,  and  the  eigenvalue  A2  =  ~7n-i  -  i^n-l  nas  an 
associated  eigenvector  x2  :=  [2"    ,2"   J 

Note  that  since  the  eigenvalues  of  G  given  by  (4.2)  are  A  =  ±1, 
independent  of  -1  <  7  <  1,  deflation  takes  place  numerous  times  during  the 
computat  ions . 


We  turn  to  the  computation  of  the  Householder  transformation  (1 .4) . 
In  oder  to  avoid  cancellation  of  significant  digits,  we  compute  {ws?ws-»-i} 
given  by  (1.6)  as  follows.    If  7S  >  0,  then  we  use  (1.6a)  and  replace 
(1.6b)  by 


ws  +  i  :-  -  <7s2    (1  +  7S) 


(1.6b') 


In    case    7S    <    0 ,     we    use     (1.6b)     and    replace     (1.6a)     by 


,,    ._  „   0-1/2 M  _  >.-l/2 

u>s  .-  asZ         (1     -    7S; 


(1.6a') 


Due  to  H  having  real-valued  elements,  the  eigenvalues  and  eigenvectors 
of  H  occur  in  complex  conjugate  pairs.   Therefore  only  zeros  of  $(#)  for 
0  <  6    <    ir    have  to  be  computed.   Moreover, 

n     _     ,9-    6, 
•(0  =  .£  Kjl2  cot(^-) 
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can  be  simplified.   Assume  that  0  <  #k  <  n    for  some  k  <  n  and  let 
^k  +  1  -  27r  -  0k  .   Then  |  (k  |  =  Kk+1 I '  and  we  obtain 


'Gi.  i_i  -9 


9^-9< 


'0^+0 


|Ck|2cot(V)    +     Kk  +  1!2cot(^^)    =     ,Ck|2(cot(V)  cot(^)) 


=  I  Ck  I 


2  sin  9 


COS  V     -     cos 


Kkl2 


sin  9 


■>»m*H^) 


(4.5) 


We  use  the  right  hand  side  of  (4.5)  in  the  computations 
we  need  to  evaluate  cot(-l)  as  well  as  cotf  n~    J  =  tan  S. 
The  contribution  from  (4.5)  to  $'(0)     is 

1     -     cos    9    cos    9U 


2lCkl2    dflVcos    01-    cos    9 J    ~    2|Ckl 


(cos 


cos  9k)2 


If  9k    =    0  then 


(4.6) 


The  stable  evaluation  of  the  right  hand  side  of  (4.6)  can  be  accomplished 
as  described  in  Table  4.1. 


Cond  i t  i  ons 


Eval uate 


cck  <  0 


cck  >  0  and  s  cs_  >  0 


cck  >  0  and  s  cs_  <  0 


1  /  -2     -2n 
J(sjf  +  s-z) 


4^+ 


2s  i  s- 


Us+s-J 

Ck    +  (      S*   V 
2s,  s-    \2s,  s-/ 


Table  4.1:   Stable  evaluation  of  ( 1 -cck ) / (c-ck )  ,  where 


c  : =  cos  9 ,  ck  : 

:=  sintV)- 


r9„+9 


cos  Pk  ,  s  :  =  sin  V  ,  sk  :=  sin  Pk  ,  s.  :=    si 


n(-V), 


The  interlacing  of  the  zeros  of  <^(A)  with  the  {Ak}j)  ,  implies  that  it 
easily  can  be  determined  whether  0    -    0  or  0    =  n    are  zeros  of  <I>(0)  .   Let 
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the  0k  ,  1  <  k  <  n  ,  be  ordered  so  that  0  <  6X     <    02     <     ...  <  0p  <  n    <  0p  +  1 
<  ...  <  0n    <    2n  .   Since  the  Ak  =  exp(i#k)  appear  in  complex  conjugate 
pairs,  we  obtain 

c6l    >   0  =>      $(0)  =   0  , 
\ep    <    n       =>   $(tt)  =  0  . 

Finally,  we  consider  the  computation  of  eigenvectors  v»  defined 
by(2.20) .   Let 

W2=i   [v«  w«  .  .  .  ,«£2]  ,      wf  6C"-  , 

Ax    =:     diag[exp(i«j  ')  ,exp(i^    )  ,  .  .  .  ,exp(iffg  ')]  ,  0    <    0-    '    <    2w     , 

(2)  (2)  (2)  (2) 

A2    =:     diag[exp(iflj  ;)  ,exp(i^  ;)  ,  .  .  .  ,exp(  i0},.£)]  »  0    <    0J   ;    <    2*     , 

«1     es^s    =!      LCi     ><2     '  •  •  •  »Cs    J        > 

w  H.  ,,         -•     iy(2>   ^(2)  /2)lT 

w2    elu's  +  i    — •      Lsj     » ^2     '  '  •  '  'sn-sj         ' 


and    A    =:     exp(i0),     0    <    9    <    2n  .       Then 

W1(I-A1HA)-1W1Hescs    =    £      (1    -    expCiC*-^)))-1    Ci^M^ 

j=l  J  J        J 

E  (1     -    expCKfl-^)))-^ 

(1)  J  J        J 

+  E  [(l-expCi^-^)))"^}1^    +     (l-explilJ))))^] 

n\  J  JJ  J  JJ 

O<0j   }<n 


28 


=  1    £  (l  +  lcot^xj^wj"  +  i  Z      (l-itanCl))^1^1' 


<f=0 


0-     =n 
J 


0<*(1)<*  9    .  ,'  V  •  ,'  '"'X 

j  2  sin(-L__jsin^^__j 


(4.7) 


.(1)  (1). 


-  i  sin  0    £   (sin(-L-)sin(-V-))^  ReCC^wj1') 

o<0:  '<* 


We  may  assume  that  close  eigenvalues  have  been  eliminated  from  A    and  A 

bv  deflation  ,  and  that  therefore  the  0-        and  0-        are  distinct.   Hence,  the 
*  J         J 

sums  over  0-        =0  and  0-        —    n    contain  at  most  one  term  each. 
Analogously  to  (4.7)  we  obtain 
-1„  H.           n^S 


W2(A2-IA)"1W2"elWs+1    £  (exp(i^2))-exp(i0))-1c|2)wj2) 

j=l        J  J   J 


1   E   (l  +  icot(|)Kj2)w^  +  1   £   (-l-fitan(|))cj2)wj2) 


<j2,=° 


,0+0 


(2) 


,(2) 


+  2 


£ 


S  1 


o<*<2><» 

J      si 


^■si^W>»f>) 


,6-e.  e+0'. 

n{-J-)       Sin(-T-) 


2e 


,(2) 


,0+0. 


(2) 


,(2) 


(2)  (2), 


E        sin  ^(sin(^_)sin(_^_))-i   i.((WwW) 


0<^<, 


+  A  sin 


,(2) 


,0+0. 


(2) 


,(2) 


(2)  (2)' 


£     cos  ^(sin(-^-)sin(-^-))-1  Re(CjWwj^) 


0<f)<, 


(4.8-) 
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The  evaluation  of  6(A)  defined  by  (2.19)  can  also  be  simplified.   We  have 


*(A)     =       E 


i<jV 


,(i) 


n-s 

+  £ 


c,(2V 


1/2 


(2). 


1     |A-exp(i0J  0|2        J=l     |A-exp(i0J  ')  | 


where ,     e.g., 

E 

j=l  |A-exp(i«j'))|2 

'   3      ? 
,f=0 

.(1) 


=  4  e  icjV/sin'cg)  +  4  e  id 


(i) 


|2/cos2(^) 


0(1)- 


+ 


(1) 


1      r     ic  ;i 


[1 

o<er'<n 


,(1) 


0+0: 


(1; 


(sin(-^-))-    +    (sin^)) 


(4.9) 


The  simplifications  of  this  section  for  the  orthogonal  eigenproblem  have 
been  implemented  in  a  Pascal  program.   Several  other  mathematically 
equivalent  forms  of  (4.7)-(4.9)  could  also  be  used.   We  have  tried  to  find 
formulas  that  avoid  unnecessary  loss  of  significant  digits. 
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5.   Numerical  examples 


We  report,  results  of  some  computed  examples  with  an  experimental 
program  for  the  orthogonal  e igenproblem .   The  program  is  written  in  Turbo 
Pascal  4.0  and  was  run  on  an  IBM  PC  AT  computer  with  unit  roundoff 
u  =  2    ~  2-10"   .   Our  code  implements  the  formulas  of  Section  4. 
Generally  very  accurate  answers  are  obtained.   Lemma  3.2  indicates, 
however,  that  a  zero  0    of  $(#)  may  be  very  close  to  a  singular  point  0     of 
$(#)  and  by  Lemma  3.3  the  difference  0-0-    has  to  be  computed  to  high 
re  1  at  i ve  accuracy  in  order  to  yield  nearly  orthogonal  eigenvectors. 
Example  5.2  below  shows  that,  indeed,  6-6-    can  be  extremely  tiny  and  that 
loss  of  accuracy  in  both  eigenvectors  and  eigenvalues  may  result.   This 
loss  of  accuracy  could  be  reduced,  e.g. ,  by  representing  6    and  0-     in  higher 
precision  arithmetic. 

In  this  section  A  G  CnXx  denotes  the  diagonal  matrix  with  the  computed 
eigenvalues  of  H  €  RnXn  as  entries,  and  W  6  CnXn  is  the  matrix  with  the 
computed  eigenvectors.   We  evaluate  the  residual  errors  HHW-WAHqq  and 
|| W  W-IHoo,  where     Hoo  denotes  the  uniform  matrix  norm. 

Example  5.1.   This  example  discusses  the  application  of  the  unitary  and 
orthogonal  e i genprob 1  ems  to  the  construction  of  Gauss-Szego  quadrature 
rules.   Consider  the  inner  product  on  the  unit  circle 


<f,g>  =  f(A)  g(A)  da(A)  , 

J|A|  =  1 


(5.1) 


vith  a  positive  measure  da(A)  .   Let  V-'k  »  0  <  k  <  n,  be  monic  orthogonal 
polynomials  with  respect  to  (5.1).   They  satisfy  a  recurrence  relation 


3] 


<A0(A)   :=  1 

V-k(A)   :=  AVk.i(A)  +  7k0k.i(A),  1  <  k  <  n 


(5.1a) 
(5.2b) 


for  some  parameters  7k  G  C  such  that  |  7k  |  <  1  for  1  <  k  <  n.   Here 
V'k.1(A)  :=  A"  Vk-i(l/^)  1S    the  "reversed  polynomial."   Let  7n  €  C  be  an 
arbitrary  complex  number  of  unit  magnitude,  and  define  ipn    by  (5.2b)  with 
k  :=  n.   Writing  the  recursions  (5.2)  (for  1  <  k  <  n)  in  matrix  form 
yields  the  unitary  matrix 


H  —  GXG2  •  •  .Gn  jGn, 


(5.3) 


whose  eigenvalues  {Ak}P,  are  the  zeros  of  V'n  •   Here  Gk  is  defined  by  7k 
according  to  (1.2)  for  1  <  k  <  n.   Hence,  the  parameters  {7k}n  are  the 

U 

Schur  parmaeters  for  H.   Let  H  =  WAW   be  a  spectral  resolution,  and  define 
the  weights  pk  :=  |e,Wek|   for  1  <  k  <  n.   Then 

f(A)da(A)  =    £     pkf(Ak)  +  en(f) 
|A|=1  k=l 


is  a  Gauss-Szego  quadrature  rule  with  respect  to  the  measure  da(A) , 
because  the  error  cn(f)  vanishes  when  f  is  any  trigonometric  polynomial  of 
degree  less  than  n.   See  [Gr2]  for  details.   The  computed  examples 
illustrate  the  case  when  all  Schur  parameters  7k  are  real  valued  and  H 
therefore  is  real  orthogonal. 

A  particularly  simple  example  is  7k  :=  0,  1  <  k  <  n,  and  7n  :=  -1. 
Then  i<k(\)     =    \k  ,     0    <    k    <    n  ,     and  t/'n(A)  =  An  -  1  ,  and  therefore 


Ak  =  exp(27ri  (k-1  )/n)   , 
.pk    =    1/n  . 


1  <  k  <  n 


(5.4) 
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n 

||HW-WA||oo 

l|WMW-I||o0 

4 

4.6-10"12 

i  .5-icr11 

8 

6.4-10"12 

3.010"11 

16 

1  ,4-lCf11 

5.4-10  n 

32 

2.910"11 

1  .8-10"10 

64 

3.9-10"11 

3.4-10"10 

These  Schur  parameters  have  been  used  for  Table  5.1.    In  the  table  "# 
defl.   close  e.v."  stands  for  number  of  deflations  due  to  close 
eigenvalues,  and  "#  defl.  small  I  C|<  I "  *s  short  for  number  of  deflations 
due  to  components  £k  of  z  of  small  magnitude.   Two  eigenvalues  are 
considered  close  if  (2.26)  is  satisfied  for  e2  :=  1-10"  ,  and  |  (k  |   is 
regarded  small  if  (2.24)  is  valid  for  cl     :—    1-10"  .   These  values  of  (^    and 
(2  are  used  in  all  computed  examples  of  this  section. 


#  defl.  close  e.v.      #  defl.  small   | (k  | 

2  0 

8  0 

24  0 

64  0 

160  0 


Table  5.1:   7k  :=  0,  1  <  k  <  n;  yn     :=    -1 


For  7k  :=  0,  1  <  k  <  n,  and  7n  :=  1,  we  obtain  the  polynomials  Vk(^)  = 
A  ,  0  <  k  <  n,  and  ^'n(^)  :  =  An  +  1.   Hence,  the  eigenvalues  are  Ak  = 
exp  (  i  7T  (2k- 1 ) /n )  ,  1  <  k  <  n,  and  the  Gauss-Szego  weights  pk  are  the  same  as 
in  (5.4).   Table  5.2  shows  computations  for  the  present  Schur  parameters, 

and  differs  from  Table  5.1  mainly  in  that  fewer  deflations  take  place. 

#  defl  .  smal  1   |  <k  | 
0 
0 

o 
o 

0 
Table  5.2:   7k  :=  0,  1  <  k  <  n;  -)n  :=  1 
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n 

HHW-WAHoo 

HWHW-I||co 

# 

defl  . 

c 1 ose    e . v 

4 

7.8-10"12 

1  .6-10*11 

0 

8 

1  .7-10"11 

4.2-10"11 

■2 

16 

3.  1-10"11 

1  .610"10 

10 

32 

4.  1-10"11 

3.5-10"10 

34 

64 

5.610"11 

7.5-10"10 

f)S 

In  Tables  5.3-5.4  we  have  chosen  7k  :=  0.8,  1  <  k  <  n.   This  makes  the 
Ak  gather  in  the  left  half  plane.   For  the  examples  of  Table  5.3  we  have 

max  Re  Xl,    <    -\.        For  the  examples  of  Table  5.4  we  obtain  max  Re  Ak  <  -  -j  . 

Ak^l       k      4  4 


n 

HHW-WAHoo 

l|WHW-I||00 

# 

def  1 

.     c 1 ose    e . v 

4 

2.7-10  12 

1  .610"11 

2 

8 

5.5-10"11 

1  .8-10"10 

8 

16 
32 

5.210"11 
3.2-10"8 

3.2-10"10 
9.3-10"8 

24 

63 

64 

3.2-10"8 

1  .610"7 

157 

#  defl 


smal  1   |  Ck 
0 

0 
0 

1 
3 


Table  5.3:   7k  :  =  0.8,  1  <  k  <  n;  7n  :=  -1 


#  defl 


n 

HHW-WAHoo 

l|WHW-I||oo 

# 

d 

ef  1 

c 1 ose    e . v 

4 

4.8-10"11 

1  .7-10"10 

0 

8 

9.4-10"11 

5.5-10-10 

2 

16 

4.8-10"10 

6.6-10"9 

10 

32 

6.310"10 

2.310"8 

34 

64 

4.210"8 

1  .9-10"7 

97 

smal  1  |  Ck  I 
0 
0 
0 
0 
1 


Table  5.4:   7k:=0.8,l<k<n;7n:=l 

In  the  last  computed  quadrature  rules  of  this  example  we  let  the  7k  ,  1 
<  k  <  n,  be  uniformly  distributed  in  the  open  interval  ]-l,l[,  and  let  7n 
be  -1  or  1  with  probability  i  each.   The  7k  are  determined  with  the  random 
number  generator  of  Pascal.   Table  5.5  shows  the  result  of  30 
e igenprobl ems  so  generated.   The  maximum,  average  and  minimum  in  Table  5.5 
are  over  all  30  e igenproblems . 
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max 

average 
mi  n 


HHW-WAHoo 
7.2-107 
5.8-10"8 
2.910"9 


||WHW-  I||oo  #  defl.  close  e.v 

2.5-10'6  30 

2.3107  26.5 

1.5-10"8  22 


#  defl  .  smal  1   |  (k 
0 
0 
0 


Table  5.5:   Uniformly  distributed  7k  e]-l,l[,  1  <  k  <  n;  uniformly 
distributed  7n  £  {-1,1}.   Max,  average  and  min  are  over 
30  e igenprobl ems  with  n  :=  32 


The  numerical  experiments  of  Table  5.5  indicate  that  for  many  choices 
of  Schur  parameters  7k  ,  the  magnitudes  | (k  |  are  not  sufficiently  small  to 
give  rise  to  frequent  deflations.   This  behavior  has  also  been  observed  in 
many  other  computed  experiments.    In  contrast,  massive  deflation  in  DC 
methods  for  symmetric  tridiagonal  matrices  often  is  caused  by  small 
components  of  the  vector  correspnding  to  z  =  [Cj<]  k_i  •        E 

Example  5.2.   This  example  suggests  that  it  might  not  be  possible  to 

increase  the  small  lower  bound  for  min  |#-#j|  of  Lemma  3.2  significantly. 

J 

The  Schur  parameters  for  Table  5.6  are  obtained  by  reversing  the  sign  of 
the  7k  ,  1  <  k  <  n ,  of  Table  5.4. 


min 


n 

l<k<n 

4 

6.6-10 

8 

8.  1-10 

16 

0* 

32 

0* 

64 

0* 

-2 


Table  5.6 


#    defl . 

# 

defl  . 

HW-WA||oo 

||WHW-I||oo 

c  1  ose 

e.v. 

sma 

ii    Kk 

6.9-10"11 

3.  1-10"10 

0 

0 

7.2-10"10 

2.5-10"8 

2 

0 

1  .2-10"7 

3.1-10"8 

10 

0 

7.210"7 

1  .9-10"6 

34 

2 

7.210"7 

2.610"6 

97 

5 

7k  :=  -0.8,  1  <  k  <  n;  7n  :=  1.   *The  matrix  has 
numerically  the  eigenvalue  X    =    1  of  multiplicity  two 
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Because  ak    =  0.6  >  0,  1  <  k  <  n,  "the  matrix  H  has  distinct  eigenvalues 
mathematically.   Numerically  two  eigenvalues  are  so  close  that  they  are 
not  distinguished  with  our  present  choice  of  e2  =  1-10"  .   A  smaller  value 
of  £2 »  such  as  €2  =  1*10"  ,  gave  in  some  numerical  experiments  larger 
residual  errors  ||HW-WA||oo  or  ||WHW-I||00.      □ 
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